setwd("/Users/naimagreen/Documents/Quant Files/GOV 2001/Replication Paper")


#replication code for gov 2001 final paper

require(MASS)
require(MatchIt)
require(stargazer)

load("exp1_rep.RData")
#load data used in exp1 and 2 as exp1 and exp2
exp1<-data.frame(x)
load("exp2_rep.Rdata")
exp2<-data.frame(x)

#show that covariates balance is not perfect
#exp1, rebuta as the treatment
exp1match<-exp1[-c(11,12)]
m.out.exp1ra <- matchit(rebuta ~life+news+politicalinterest+age+education+income+partymember+prowest+male, data = exp1match, method = "genetic")
summary(m.out.exp1ra)
t<-summary(m.out.exp1ra)
stargazer(t$sum.all)
######################################################################################
#replicate table 1,2,4,5 for appendix
group1<-subset(exp1,group==1)
group2<-subset(exp1,group==2)
group3<-subset(exp1, group==3)
group4<-subset(exp1,group==4)
group5<-subset(exp1,group==5)
#table1
t.test.se <- function(var1,var2)
{
  m1<-mean(var1)
  m2<-mean(var2)
  s1<-sd(var1)
  s2<-sd(var2)
  n1<-length(var1)
  n2<-length(var2)
  se <- sqrt( (s1^2/n1) + (s2^2/n2) )
  meandiff<-m1-m2
  df <- ( (s1^2/n1 + s2^2/n2)^2 )/( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )
  t <- (m1-m2)/se
  result<-c(meandiff,se,2*pt(-abs(t),df))
  names(result)<-c("difference in mean", "standard error","p-value")
  return(result)
}
t.test.se(group1$safetytrust,group2$safetytrust)
t.test.se(group1$equalitytrust,group2$equalitytrust)
t.test.se(group1$polisystem,group2$polisystem)
#table1,second part
t.test.se(group3$rumorabelief,group2$rumorabelief)
t.test.se(group3$rumorbbelief,group2$rumorbbelief)
t.test.se(group3$safetytrust,group2$safetytrust)
t.test.se(group3$equalitytrust,group2$equalitytrust)
t.test.se(group3$polisystem,group2$polisystem)
t.test.se(group4$rumorabelief,group2$rumorabelief)
t.test.se(group4$rumorbbelief,group2$rumorbbelief)
t.test.se(group4$safetytrust,group2$safetytrust)
t.test.se(group4$equalitytrust,group2$equalitytrust)
t.test.se(group4$polisystem,group2$polisystem)
t.test.se(group5$rumorabelief,group2$rumorabelief)
t.test.se(group5$rumorbbelief,group2$rumorbbelief)
t.test.se(group5$safetytrust,group2$safetytrust)
t.test.se(group5$equalitytrust,group2$equalitytrust)
t.test.se(group5$polisystem,group2$polisystem)
#table1,third part
t.test.se(group3$safetytrust,group1$safetytrust)
t.test.se(group3$equalitytrust,group1$equalitytrust)
t.test.se(group3$polisystem,group1$polisystem)
t.test.se(group4$safetytrust,group1$safetytrust)
t.test.se(group4$equalitytrust,group1$equalitytrust)
t.test.se(group4$polisystem,group1$polisystem)
t.test.se(group5$safetytrust,group1$safetytrust)
t.test.se(group5$equalitytrust,group1$equalitytrust)
t.test.se(group5$polisystem,group1$polisystem)


########################################################################
#replication:table 3 in original model, and table 3 in logit model
ps1 <- polr(as.factor(safetytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="probit")
ps1logit <- polr(as.factor(safetytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="logistic")
et1 <- polr(as.factor(equalitytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="probit")
et1logit <- polr(as.factor(equalitytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="logistic")
psys1 <- polr(as.factor(polisystem) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="probit")
psyslogit <- polr(as.factor(polisystem) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp1, Hess = TRUE, method ="logistic")
stargazer(ps1,ps1logit,et1,et1logit,psys1,psyslogit)
##############################################################################################
#replication:table 6 in original model, and table 6 in logit model
cp1 <- polr(as.factor(protectiontrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="probit")
cp1logit <- polr(as.factor(protectiontrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="logistic")
psys2 <- polr(as.factor(polisystem) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="probit")
psys2logit <- polr(as.factor(polisystem) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="logistic")
of1 <- polr(as.factor(officialstrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="probit")
of1logit <- polr(as.factor(officialstrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = exp2, Hess = TRUE, method ="logistic")
stargazer(cp1,cp1logit,of1,of1logit,psys2,psyslogit)
####################################################################################

# Proportional Odds Plot for Equal Treatment Trust

load("exp1_rep.Rdata")

require(ordinal)
require(foreign)
require(ggplot2)
require(Hmisc)
require(reshape2)


## fit ordered progit model and store results
equaltreatment <- polr(as.factor(equalitytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(equaltreatment)




#plot 

## store table
(ctable7 <- coef(summary(equaltreatment)))

## calculate and store p values
p7 <- pnorm(abs(ctable7[, "t value"]), lower.tail = FALSE) * 2

## combined table
(ctable7 <- cbind(ctable7, "p value" = p7))

(ci7 <- confint(equaltreatment)) # default method gives profiled CIs

confint.default(equaltreatment) # CIs assuming normality

## odds ratios
exp(coef(equaltreatment))

## OR and CI
exp(cbind(OR = coef(equaltreatment), ci7))

sf <- function(y) {
  c('Y>=1' = qlogis(mean(y >= 1)),
    'Y>=2' = qlogis(mean(y >= 2)),
    'Y>=3' = qlogis(mean(y >= 3)),
    'Y>=4' = qlogis(mean(y >= 4)),
    'Y>=5' = qlogis(mean(y >= 5)),
    'Y>=6' = qlogis(mean(y >= 6)),
    'Y>=7' = qlogis(mean(y >= 7)))
}


(s7 <- with(x, summary(as.numeric(equalitytrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, fun=sf)))

## Setting difference based on normalizing difference from Y = 1 to Y = 2

s7[, 7] <- s7[, 7] - s7[, 3]
s7[, 6] <- s7[, 6] - s7[, 3]
s7[, 5] <- s7[, 5] - s7[, 3]
s7[, 4] <- s7[, 4] - s7[, 3]
s7[, 3] <- s7[, 3] - s7[, 3]
s7 # print

#Final Plot

pdf("huangPOA8.pdf")
plot(s7, which=1:7, pch=c(1,7,8,11,17,23), bg = "red", cex = .7,  xlab='probit', main='Testing the Proportional Odds Assumption \n Trust in Equal Treatment', xlim=c(-5,0))
dev.off()

##################################################################

#generate the histogram of levels of education to show violation of POA
exp1edu1<-subset(exp1,education==1)
exp1edu2<-subset(exp1,education==2)
exp1edu3<-subset(exp1,education==3)
exp1edu4<-subset(exp1,education==4)
exp1edu5<-subset(exp1,education==5)
exp1edu6<-subset(exp1,education==6)
#equalitytrust proportion for each education level
etedu1<-as.data.frame(table(exp1edu1$equalitytrust))
etedu2<-as.data.frame(table(exp1edu2$equalitytrust))
etedu3<-as.data.frame(table(exp1edu3$equalitytrust))
etedu4<-as.data.frame(table(exp1edu4$equalitytrust))
etedu5<-as.data.frame(table(exp1edu5$equalitytrust))
etedu6<-as.data.frame(table(exp1edu6$equalitytrust))
etedu2$percent<- prop.table(etedu2$Freq)
etedu2
etedu3$percent<- prop.table(etedu3$Freq)
etedu3
etedu4$percent<- prop.table(etedu4$Freq)
etedu4
etedu5$percent<- prop.table(etedu5$Freq)
etedu5
etedu6$percent<- prop.table(etedu3$Freq)
etedu6
#polisystem for each edu level, exp1
ps1edu1<-as.data.frame(table(exp1edu1$polisystem))
ps1edu2<-as.data.frame(table(exp1edu2$polisystem))
ps1edu3<-as.data.frame(table(exp1edu3$polisystem))
ps1edu4<-as.data.frame(table(exp1edu4$polisystem))
ps1edu5<-as.data.frame(table(exp1edu5$polisystem))
ps1edu6<-as.data.frame(table(exp1edu6$polisystem))
ps1edu1$percent<- prop.table(ps1edu1$Freq)
ps1edu1
ps1edu2$percent<- prop.table(ps1edu2$Freq)
ps1edu2
ps1edu3$percent<- prop.table(ps1edu3$Freq)
ps1edu3
ps1edu4$percent<- prop.table(ps1edu4$Freq)
ps1edu4
ps1edu5$percent<- prop.table(ps1edu5$Freq)
ps1edu5
ps1edu6$percent<- prop.table(ps1edu6$Freq)
ps1edu6
#subset exp2 data based on edu level
exp2edu1<-subset(exp2,education==1)
exp2edu2<-subset(exp2,education==2)
exp2edu3<-subset(exp2,education==3)
exp2edu4<-subset(exp2,education==4)
exp2edu5<-subset(exp2,education==5)
exp2edu6<-subset(exp2,education==6)
#citizenprotection proportion for each education level
cpedu1<-as.data.frame(table(exp2edu1$protectiontrust))
cpedu2<-as.data.frame(table(exp2edu2$protectiontrust))
cpedu3<-as.data.frame(table(exp2edu3$protectiontrust))
cpedu4<-as.data.frame(table(exp2edu4$protectiontrust))
cpedu5<-as.data.frame(table(exp2edu5$protectiontrust))
cpedu6<-as.data.frame(table(exp2edu6$protectiontrust))
cpedu2$percent<- prop.table(cpedu2$Freq)
cpedu2
cpedu3$percent<- prop.table(cpedu3$Freq)
cpedu3
cpedu4$percent<- prop.table(cpedu4$Freq)
cpedu4
cpedu5$percent<- prop.table(cpedu5$Freq)
cpedu5
cpedu6$percent<- prop.table(cpedu6$Freq)
cpedu6
#polisystem for each edu level, exp2
ps2edu1<-as.data.frame(table(exp2edu1$polisystem))
ps2edu2<-as.data.frame(table(exp2edu2$polisystem))
ps2edu3<-as.data.frame(table(exp2edu3$polisystem))
ps2edu4<-as.data.frame(table(exp2edu4$polisystem))
ps2edu5<-as.data.frame(table(exp2edu5$polisystem))
ps2edu6<-as.data.frame(table(exp2edu6$polisystem))
ps2edu1$percent<- prop.table(ps2edu1$Freq)
ps2edu1
ps2edu2$percent<- prop.table(ps2edu2$Freq)
ps2edu2
ps2edu3$percent<- prop.table(ps2edu3$Freq)
ps2edu3
ps2edu4$percent<- prop.table(ps2edu4$Freq)
ps2edu4
ps2edu5$percent<- prop.table(ps2edu5$Freq)
ps2edu5
ps2edu6$percent<- prop.table(ps2edu6$Freq)
ps1edu6


########################################################################################
########################################################################################
########################################################################################
# Full Replication of Huang's Results
#
########################################################################################
########################################################################################
########################################################################################

## Experiment 1 Tables and Plots
load("exp1_rep.Rdata")

#table2
rumor<-rbind(group2,group3,group4,group5)
rumor[sapply(rumor, is.factor)] <- lapply(rumor[sapply(rumor, is.factor)], as.numeric)
rumor$rumorabelief<-as.factor(rumor$rumorabelief)
rumora<-polr(rumorabelief~ rebuta + rebutb + rebut2 + news+politicalinterest + life + prowest +age+male+education+income+partymember, data= rumor, method="probit")
rumor$rumorbbelief<-as.factor(rumor$rumorbbelief)
rumorb<-polr(rumorbbelief~ rebuta + rebutb + rebut2 + news+politicalinterest + life + prowest +age+male+education+income+partymember, data= rumor, method="probit")
stargazer(rumora,rumorb)
#====================================================================
load("exp2_rep.Rdata")
## difference in means from Table 4
#First, some brainstorming on how to find the differences in means
#row 1, column 3
protect_RC_test <- t.test(x$protectiontrust[x$group == 2], x$protectiontrust[x$group == 1], var.equal = F, alternative = "less")
protect_RC_diff <-as.data.frame(protect_RC_test[5])[1,1]-as.data.frame(protect_RC_test[5])[2,1]
protect_RC_diff
protect_RC_SE <- protect_RC_diff/(as.data.frame(protect_RC_test[1])[1,1])
protect_RC_SE
protect_RC_pval <- as.data.frame(protect_RC_test[3])[1,1]
protect_RC_pval
#row 1, column 4
fam_RC_test <- t.test(x$officialstrust[x$group == 2], x$officialstrust[x$group == 1], var.equal = F, alternative = "less")
fam_RC_diff <-as.data.frame(fam_RC_test[5])[1,1]-as.data.frame(fam_RC_test[5])[2,1]
fam_RC_diff
fam_RC_SE <- fam_RC_diff/(as.data.frame(fam_RC_test[1])[1,1])
fam_RC_SE
fam_RC_pval <- as.data.frame(fam_RC_test[3])[1,1]
fam_RC_pval
# Now, to make Table 4
#First, a matrix of the various groups to be used
left <- c(2, 3, 4, 5, 3, 4, 5)
right <- c(1, 2, 2, 2, 1, 1, 1)
order <- as.matrix(cbind(left, right))
order
# Now, finding the values for each column
pass_belief <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
  }
  else { pass_belief[i] <- NA}
}
pass_belief
ferrari_belief <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
  }
  else { ferrari_belief[i] <- NA}
}
ferrari_belief
protect_trust <- rep(NA, 7)
for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
}
protect_trust
fam_trust <- rep(NA, 7)
for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
}
fam_trust
system_trust <- rep(NA, 7)
for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
}
system_trust
# Adding a column of titles
titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")
# Combining all columns together to make the table
table4_diff <- cbind(titles, pass_belief, ferrari_belief, protect_trust, fam_trust, system_trust)
table4_diff
# Doing the same for Standard Errors
pass_belief_SE <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
    pass_belief_SE[i] <- pass_belief[i]/(as.data.frame(pass_belief_test[1])[1,1])
  }
  else { pass_belief_SE[i] <- NA}
}
pass_belief_SE
ferrari_belief_SE <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
    ferrari_belief_SE[i] <- ferrari_belief[i]/(as.data.frame(ferrari_belief_test[1])[1,1])
  }
  else { ferrari_belief_SE[i] <- NA}
}
ferrari_belief_SE
protect_trust_SE <- rep(NA, 7)
for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
  protect_trust_SE[i] <-  protect_trust[i]/(as.data.frame(protect_trust_test[1])[1,1])
}
protect_trust_SE
fam_trust_SE <- rep(NA, 7)
for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
  fam_trust_SE[i] <-  fam_trust[i]/(as.data.frame(fam_trust_test[1])[1,1])
}
fam_trust_SE
system_trust_SE <- rep(NA, 7)
for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
  system_trust_SE[i] <-  system_trust[i]/(as.data.frame(system_trust_test[1])[1,1])
}
system_trust_SE
# Adding a column of titles
titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")
# Combining all columns together to make the table
table4_SE<- cbind(titles, pass_belief_SE, ferrari_belief_SE, protect_trust_SE, fam_trust_SE, system_trust_SE)
table4_SE
#Finally, let's make a table for statistical significance
pass_belief_SS <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
    pass_belief_pval <- as.data.frame(pass_belief_test[3])[1,1]
    if (pass_belief_pval <= .01) {pass_belief_SS[i] <- "***"}
    else if (pass_belief_pval > .01 & pass_belief_pval <= .05){pass_belief_SS[i] <- "**"}
    else if (pass_belief_pval > .05 & pass_belief_pval <= .1){pass_belief_SS[i] <- "*"}
    else {pass_belief_SS[i] <- "not SS"}
  }
  else { pass_belief_SE[i] <- NA}
}
pass_belief_SS
ferrari_belief_SS <- rep(NA, 7)
for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
    ferrari_belief_pval <- as.data.frame(ferrari_belief_test[3])[1,1]
    if (ferrari_belief_pval <= .01) {ferrari_belief_SS[i] <- "***"}
    else if (ferrari_belief_pval > .01 & ferrari_belief_pval <= .05){ferrari_belief_SS[i] <- "**"}
    else if (ferrari_belief_pval > .05 & ferrari_belief_pval <= .1){ferrari_belief_SS[i] <- "*"}
    else {ferrari_belief_SS[i] <- "not SS"}
  }
  else { ferrari_belief_SE[i] <- NA}
}
ferrari_belief_SS
protect_trust_SS <- rep(NA, 7)
protect_trust_pval <- rep(NA, 7)
for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
  protect_trust_pval[i] <- as.data.frame(protect_trust_test[3])[1,1]
  if (protect_trust_pval[i] <= .01) {protect_trust_SS[i] <- "***"}
  else if (protect_trust_pval[i] > .01 & protect_trust_pval[i] <= .05){protect_trust_SS[i] <- "**"}
  else if (protect_trust_pval[i] > .05 & protect_trust_pval[i] <= .1){protect_trust_SS[i] <- "*"}
  else {protect_trust_SS[i] <- "not SS"}}
protect_trust_SS
protect_trust_pval
fam_trust_SS <- rep(NA, 7)
fam_trust_pval <- rep(NA, 7)
for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
  fam_trust_pval[i] <- as.data.frame(fam_trust_test[3])[1,1]
  if (fam_trust_pval[i] <= .01) {fam_trust_SS[i] <- "***"}
  else if (fam_trust_pval[i] > .01 & fam_trust_pval[i] <= .05){fam_trust_SS[i] <- "**"}
  else if (fam_trust_pval[i] > .05 & fam_trust_pval[i] <= .1){fam_trust_SS[i] <- "*"}
  else {fam_trust_SS[i] <- "not SS"}}
fam_trust_SS
fam_trust_pval
system_trust_SS <- rep(NA, 7)
system_trust_pval <- rep(NA, 7)
for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
  system_trust_pval[i] <- as.data.frame(system_trust_test[3])[1,1]
  if (system_trust_pval[i] <= .01) {protect_trust_SS[i] <- "***"}
  else if (system_trust_pval[i] > .01 & system_trust_pval[i] <= .05){system_trust_SS[i] <- "**"}
  else if (system_trust_pval[i] > .05 & system_trust_pval[i] <= .1){system_trust_SS[i] <- "*"}
  else {system_trust_SS[i] <- "not SS"}}
system_trust_SS
system_trust_pval
### aaaaaand, again
# Adding a column of titles
titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")
# Combining all columns together to make the table
table4_SS<- cbind(titles, pass_belief_SS, ferrari_belief_SS, protect_trust_SS, fam_trust_SS, system_trust_SS)
table4_SS
#### this table has problemssssssss!
## Table 5
# part 1
## fit ordered progit model and store results
table5_model1 <- polr(as.factor(rumorabelief) ~ rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table5_model1)
## fit ordered progit model and store results
table5_model2 <- polr(as.factor(rumorbbelief) ~ rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
stargazer(table5_model1,table5_model2)



## Experiment 2 Tables and Plots
load("exp2_rep.Rdata")
## Figure 3

# a

belief_table <- subset(x, select = c(group, rumorabelief, rumorbbelief))

belief_table <- na.omit(belief_table, group == 1)

belief_table

table <- belief_table[order(belief_table$group),]
table

mean(table$rumorabelief[table$group == 2])

group <- c(2, 3, 4, 5)
rumora <- c(mean(table$rumorabelief[table$group == 2]), mean(table$rumorabelief[table$group == 3]), mean(table$rumorabelief[table$group == 4]), mean(table$rumorabelief[table$group == 5]))
rumorb <- c(mean(table$rumorbbelief[table$group == 2]), mean(table$rumorbbelief[table$group == 3]), mean(table$rumorbbelief[table$group == 4]), mean(table$rumorbbelief[table$group == 5]))

df_fig3a <- data.frame(group, rumora, rumorb)
df_fig3a  
fig3a <- t(df_fig3a)[-1,]  

pdf("huang_fig3a.pdf")
barplot(fig3a, beside = T, ylim = c(3.6, 5.1), xpd = F, main = "Belief in the Rumors", ylab = "Belief")
box()
text(1.5, 5.067, "5.03")
text(2.47, 4.51, "4.46")
text(4.45, 4.4, "4.35")
text(5.46, 4.4, "4.34")
text(7.46, 4.76, "4.71")
text(8.48, 3.84, "3.78")
text(10.48, 4.22, "4.17")
text(11.5, 3.9, "3.85")
axis(1, at=c(2, 5, 8, 11), labels = c("Rumors", "Rebut passport", "Rebut Ferrari", "Rebut both"))
legend(10.55, 5.087,  cex = 0.5, c("Passport Rumor", "Ferrari Rumor"), fill = c("gray39", "lightgrey"))
dev.off()

# b

protection<- subset(x, select = c(group, protectiontrust))

protection

table2 <- protection[order(protection$group),]

table2

group2 <- c(1, 2, 3, 4, 5)
protection2 <- c(mean(table2$protectiontrust[table2$group == 1]), mean(table2$protectiontrust[table2$group == 2]), mean(table2$protectiontrust[table2$group == 3]), mean(table2$protectiontrust[table2$group == 4]), mean(table2$protectiontrust[table2$group == 5]))

df_fig3b <- data.frame(group2, protection2)
df_fig3b  

fig3b <- t(df_fig3b)[-1,]  
fig3b

png(file="huang_fig3b.png",width=800,height=600,res=100)
barplot(fig3b, beside = T, ylim = c(3.1, 4.6), col = "gray39", xpd = F, main = "Trust on Citizen Protection", ylab = "Trust")
box()
axis(1, at=c(.7,1.9,3.1,4.3,5.5), labels = c("Control", "Rumors", "Rebut passport", "Rebut Ferrari", "Rebut both"))
text(.7, 4.45, "4.41")
text(1.9, 4.17, "4.13")
text(3.1, 4.54, "4.50")
text(4.3, 4.5, "4.46")
text(5.5, 4.45, "4.41")
dev.off()

# c

families <- subset(x, select = c(group, officialstrust))

families

table3 <- families[order(families$group),]

table3

group3 <- c(1, 2, 3, 4, 5)
families2 <- c(mean(table3$officialstrust[table3$group == 1]), mean(table3$officialstrust[table3$group == 2]), mean(table3$officialstrust[table3$group == 3]), mean(table3$officialstrust[table3$group == 4]), mean(table3$officialstrust[table3$group == 5]))

df_fig3c <- data.frame(group3, families2)
df_fig3c  

fig3c <- t(df_fig3c)[-1,]  
fig3c

png(file="huang_fig3c.png",width=800,height=600,res=100)
barplot(fig3c, beside = T, ylim = c(1.5, 3.0), col = "gray39", xpd = F, main = "Trust on Officials' Families", ylab = "Trust")
box()
axis(1, at=c(.7,1.9,3.1,4.3,5.5), labels = c("Control", "Rumors", "Rebut passport", "Rebut Ferrari", "Rebut both"))
text(.7, 2.48, "2.44")
text(1.9, 2.57, "2.53")
text(3.1, 2.7, "2.66")
text(4.3, 2.9, "2.86")
text(5.5, 2.56, "2.52")
dev.off()

# d

polsystem <- subset(x, select = c(group, polisystem))

polsystem

table4 <- polsystem[order(polsystem$group),]

table4

group4 <- c(1, 2, 3, 4, 5)
polsystem2 <- c(mean(table4$polisystem[table4$group == 1]), mean(table4$polisystem[table4$group == 2]), mean(table4$polisystem[table4$group == 3]), mean(table4$polisystem[table4$group == 4]), mean(table4$polisystem[table4$group == 5]))

df_fig3d <- data.frame(group4, polsystem2)
df_fig3d  

fig3d <- t(df_fig3d)[-1,]  
fig3d

png(file="huang_fig3d.png",width=800,height=600,res=100)
barplot(fig3d, beside = T, ylim = c(3.4, 4.9), col = "gray39", xpd = F, main = "Trust on Political System", ylab = "Trust")
box()
axis(1, at=c(.7,1.9,3.1,4.3,5.5), labels = c("Control", "Rumors", "Rebut passport", "Rebut Ferrari", "Rebut both"))
text(.7, 4.79, "4.75")
text(1.9, 4.5, "4.46")
text(3.1, 4.68, "4.64")
text(4.3, 4.63, "4.59")
text(5.5, 4.88, "4.84")
dev.off()

## difference in means from Table 4

#First, some brainstorming on how to find the differences in means
#row 1, column 3

protect_RC_test <- t.test(x$protectiontrust[x$group == 2], x$protectiontrust[x$group == 1], var.equal = F, alternative = "less")
protect_RC_diff <-as.data.frame(protect_RC_test[5])[1,1]-as.data.frame(protect_RC_test[5])[2,1]
protect_RC_diff

protect_RC_SE <- protect_RC_diff/(as.data.frame(protect_RC_test[1])[1,1])
protect_RC_SE

protect_RC_pval <- as.data.frame(protect_RC_test[3])[1,1]  
protect_RC_pval

#row 1, column 4

fam_RC_test <- t.test(x$officialstrust[x$group == 2], x$officialstrust[x$group == 1], var.equal = F, alternative = "less")
fam_RC_diff <-as.data.frame(fam_RC_test[5])[1,1]-as.data.frame(fam_RC_test[5])[2,1]
fam_RC_diff

fam_RC_SE <- fam_RC_diff/(as.data.frame(fam_RC_test[1])[1,1])
fam_RC_SE

fam_RC_pval <- as.data.frame(fam_RC_test[3])[1,1]  
fam_RC_pval

# Now, to make Table 4

#First, a matrix of the various groups to be used

left <- c(2, 3, 4, 5, 3, 4, 5)
right <- c(1, 2, 2, 2, 1, 1, 1)

order <- as.matrix(cbind(left, right))
order

# Now, finding the values for each column

pass_belief <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
  }
  
  else { pass_belief[i] <- NA}
}
pass_belief

ferrari_belief <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
  }
  
  else { ferrari_belief[i] <- NA}
}
ferrari_belief

protect_trust <- rep(NA, 7)

for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
}

protect_trust

fam_trust <- rep(NA, 7)

for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
}

fam_trust

system_trust <- rep(NA, 7)

for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
}

system_trust

# Adding a column of titles

titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")

# Combining all columns together to make the table
table4_diff <- cbind(titles, pass_belief, ferrari_belief, protect_trust, fam_trust, system_trust)
table4_diff

# Doing the same for Standard Errors

pass_belief_SE <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
    pass_belief_SE[i] <- pass_belief[i]/(as.data.frame(pass_belief_test[1])[1,1])
  }
  
  else { pass_belief_SE[i] <- NA}
}
pass_belief_SE

ferrari_belief_SE <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
    ferrari_belief_SE[i] <- ferrari_belief[i]/(as.data.frame(ferrari_belief_test[1])[1,1])
  }
  
  else { ferrari_belief_SE[i] <- NA}
}
ferrari_belief_SE

protect_trust_SE <- rep(NA, 7)

for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
  protect_trust_SE[i] <-  protect_trust[i]/(as.data.frame(protect_trust_test[1])[1,1])
}

protect_trust_SE

fam_trust_SE <- rep(NA, 7)

for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
  fam_trust_SE[i] <-  fam_trust[i]/(as.data.frame(fam_trust_test[1])[1,1])
}

fam_trust_SE

system_trust_SE <- rep(NA, 7)

for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
  system_trust_SE[i] <-  system_trust[i]/(as.data.frame(system_trust_test[1])[1,1])
}

system_trust_SE

# Adding a column of titles

titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")

# Combining all columns together to make the table
table4_SE<- cbind(titles, pass_belief_SE, ferrari_belief_SE, protect_trust_SE, fam_trust_SE, system_trust_SE)
table4_SE

#Finally, let's make a table for statistical significance

pass_belief_SS <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    pass_belief_test <- t.test(x$rumorabelief[x$group == order[i,1]], x$rumorabelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    pass_belief[i] <- as.data.frame(pass_belief_test[5])[1,1] - as.data.frame(pass_belief_test[5])[2,1]
    pass_belief_pval <- as.data.frame(pass_belief_test[3])[1,1] 
    if (pass_belief_pval <= .01) {pass_belief_SS[i] <- "***"}
    else if (pass_belief_pval > .01 & pass_belief_pval <= .05){pass_belief_SS[i] <- "**"}
    else if (pass_belief_pval > .05 & pass_belief_pval <= .1){pass_belief_SS[i] <- "*"}
    else {pass_belief_SS[i] <- "not SS"}
  }
  
  else { pass_belief_SE[i] <- NA}
}
pass_belief_SS

ferrari_belief_SS <- rep(NA, 7)

for (i in 1:7){
  if (i >= 2 & i<=4){
    ferrari_belief_test <- t.test(x$rumorbbelief[x$group == order[i,1]], x$rumorbbelief[x$group == order[i,2]], var.equal = F, alternative = "less")
    ferrari_belief[i] <- as.data.frame(ferrari_belief_test[5])[1,1] - as.data.frame(ferrari_belief_test[5])[2,1]
    ferrari_belief_pval <- as.data.frame(ferrari_belief_test[3])[1,1] 
    if (ferrari_belief_pval <= .01) {ferrari_belief_SS[i] <- "***"}
    else if (ferrari_belief_pval > .01 & ferrari_belief_pval <= .05){ferrari_belief_SS[i] <- "**"}
    else if (ferrari_belief_pval > .05 & ferrari_belief_pval <= .1){ferrari_belief_SS[i] <- "*"}
    else {ferrari_belief_SS[i] <- "not SS"}
  }
  
  else { ferrari_belief_SE[i] <- NA}
}
ferrari_belief_SS

protect_trust_SS <- rep(NA, 7)
protect_trust_pval <- rep(NA, 7)

for (i in 1:7){
  protect_trust_test <- t.test(x$protectiontrust[x$group == order[i,1]], x$protectiontrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  protect_trust[i] <- as.data.frame(protect_trust_test[5])[1,1]-as.data.frame(protect_trust_test[5])[2,1]
  protect_trust_pval[i] <- as.data.frame(protect_trust_test[3])[1,1] 
  if (protect_trust_pval[i] <= .01 | protect_trust_pval[i] >= .99) {protect_trust_SS[i] <- "***"}
  else if (protect_trust_pval[i] > .01 & protect_trust_pval[i] <= .05 | protect_trust_pval[i] < .99 & protect_trust_pval[i] >= .95){protect_trust_SS[i] <- "**"}
  else if (protect_trust_pval[i] > .05 & protect_trust_pval[i] <= .1 | protect_trust_pval[i] < .95 & protect_trust_pval[i] >= .9 ){protect_trust_SS[i] <- "*"}
  else {protect_trust_SS[i] <- "not SS"}}

protect_trust_SS
protect_trust_pval

fam_trust_SS <- rep(NA, 7)
fam_trust_pval <- rep(NA, 7)

for (i in 1:7){
  fam_trust_test <- t.test(x$officialstrust[x$group == order[i,1]], x$officialstrust[x$group == order[i,2]], var.equal = F, alternative = "less")
  fam_trust[i] <- as.data.frame(fam_trust_test[5])[1,1]-as.data.frame(fam_trust_test[5])[2,1]
  fam_trust_pval[i] <- as.data.frame(fam_trust_test[3])[1,1] 
  if (fam_trust_pval[i] <= .01 | fam_trust_pval[i] >= .99) {fam_trust_SS[i] <- "***"}
  else if (fam_trust_pval[i] > .01 & fam_trust_pval[i] <= .05 | fam_trust_pval[i] < .99 & fam_trust_pval[i] >= .95){fam_trust_SS[i] <- "**"}
  else if (fam_trust_pval[i] > .05 & fam_trust_pval[i] <= .1 | fam_trust_pval[i] < .95 & fam_trust_pval[i] >= .9){fam_trust_SS[i] <- "*"}
  else {fam_trust_SS[i] <- "not SS"}}

fam_trust_SS
fam_trust_pval

system_trust_SS <- rep(NA, 7)
system_trust_pval <- rep(NA, 7)


for (i in 1:7){
  system_trust_test <- t.test(x$polisystem[x$group == order[i,1]], x$polisystem[x$group == order[i,2]], var.equal = F, alternative = "less")
  system_trust[i] <- as.data.frame(system_trust_test[5])[1,1]-as.data.frame(system_trust_test[5])[2,1]
  system_trust_pval[i] <- as.data.frame(system_trust_test[3])[1,1] 
  if (system_trust_pval[i] <= .01 | system_trust_pval[i] >= .99) {protect_trust_SS[i] <- "***"}
  else if (system_trust_pval[i] > .01 & system_trust_pval[i] <= .05 | system_trust_pval[i] < .99 & system_trust_pval[i] >= .95){system_trust_SS[i] <- "**"}
  else if (system_trust_pval[i] > .05 & system_trust_pval[i] <= .1 | system_trust_pval[i] < .95 & system_trust_pval[i] >= .9){system_trust_SS[i] <- "*"}
  else {system_trust_SS[i] <- "not SS"}}

system_trust_SS
system_trust_pval

# Adding a column of titles

titles <- c("Rumors-Control", "Rebut passport - Rumors", "Rebut Ferrari - Rumors", "Rebut both - Rumors", "Rebut passport - Control", "Rebut Ferrari - Control", "Rebut both - Control")

# Combining all columns together to make the table
table4_SS<- cbind(titles, pass_belief_SS, ferrari_belief_SS, protect_trust_SS, fam_trust_SS, system_trust_SS)
table4_SS


## Table 5

# part 1


require(MASS)
require(stargazer)

?polr

## fit ordered progit model and store results 
table5_model1 <- polr(as.factor(rumorabelief) ~ rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table5_model1)
stargazer(table5_model1)

## fit ordered progit model and store results
table5_model2 <- polr(as.factor(rumorbbelief) ~ rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table5_model2)
stargazer(table5_model2)

## fit ordered progit model and store results
table6_model1 <- polr(as.factor(protectiontrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table6_model1)
stargazer(table6_model1)

## fit ordered progit model and store results
table6_model2 <- polr(as.factor(officialstrust) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table6_model2)
stargazer(table6_model2)

## fit ordered progit model and store results
table6_model3 <- polr(as.factor(polisystem) ~ rumor + rebuta + rebutb + rebut2 + news + politicalinterest + life + prowest + age + male + education + income + partymember, data = x, Hess = TRUE, method = c("probit"))
## view a summary of the model
summary(table6_model3)
stargazer(table6_model3)

##########################################

savehistory("C:/Users/naimagreen/Desktop/replicationfinalcodes.R")